function calcUserKernel(R)

% R: N by M rating matrix, N = #users, M = #movies

N = size(R, 1);     % #users
M = size(R, 2);     % #movies

K_u = zeros(N, N);

% The generator matrix
W = zeros(N,N);
% 
% % Mean rating of each User
% U_mean = zeros(N, 1);
% for i = 1 : N
%     nzInd = (R(i,:) ~= 0);
%     U_mean(i) = mean(R(i, nzInd));
% end

% Compute the generator matrix
for i = 1 : N
    for j = i : N 
        
        %******************************************************************
        % Preparing for diffusion kernel
        
%         % Find the subset of movies that both users rated
%         comInd = find((R(i,:) ~= 0) & (R(j,:) ~= 0));
%         if (numel(comInd) == 0)
%             W(i, j) = 0;
%             continue;
%         end
%         % The Pearson r Correlation
%         numerator = 0;
%         denom_sum1 = 0;
%         denom_sum2 = 0;
%         for s = 1 : numel(comInd)
%             numerator = numerator + (R(i,comInd(s)) - U_mean(i)) * (R(j, comInd(s)) - U_mean(j));
%             denom_sum1 = denom_sum1 + (R(i,comInd(s)) - U_mean(i))^2;
%             denom_sum2 = denom_sum2 + (R(j,comInd(s)) - U_mean(j))^2;
%         end
%         denom = sqrt(denom_sum1 * denom_sum2);
%         if(denom == 0)  % FIXME: should I assign 0 in this case?
%             W(i,j) = 0;
%             continue;
%         end
%         W(i, j) = numerator/denom;

        % Cosine similarity
        W(i,j) = R(i,:) * R(j,:)'/ (norm(R(i,:)) * norm(R(j,:)));

       %*******************************************************************
  
    end
    W(:,i) = W(i,:);
end


% Construct the epsilon-neighborhood graph
thres = 0.3;
for i = 1 : N
    ind = W(i,:) >= thres;
    W(i,:) = 0;
    W(i,ind) = 1;
    W(i,i) = 0; % FIXME: necessary?
    W(:,i) = W(i,:);
end

%*******************************************************************
% Exponential kernel
K_u = expm(0.2 * W);
K_u_inv = inv(K_u);

save('u1_user_exp_kernel.mat', 'K_u');
save('u1_user_exp_kernel_inv.mat', 'K_u_inv');

%*******************************************************************
% Diffusion kernel

% beta = 0.005

% % The degree matrix
% Deg = diag(sum(W));
% 
% % Negation of graph laplacian
% L_u = W - Deg;
% 
% K_u = expm(beta * L_u);
% K_u_inv = inv(K_u);
% 
% save('u1_user_diff_kernel.mat', 'K_u');
% save('u1_user_diff_kernel_inv.mat', 'K_u_inv');

%*******************************************************************
% Commute Time kernel

% The degree matrix
% Deg = diag(sum(W));
% 
% % The graph laplacian
% L_u = Deg - W;
% 
% K_u = pinv(L_u);
% K_u_inv = inv(K_u);

% save('u1_user_ct_kernel.mat', 'K_u');
% save('u1_user_ct_kernel_inv.mat', 'K_u_inv');

%*******************************************************************
% Regularized laplacian kernel

% sigma = 5;
% 
% % The degree matrix
% Deg = diag(sum(W));
% 
% % The graph laplacian
% L_u = Deg - W;
% 
% % Regularized laplacian
% L_u_reg = Deg^(-0.5) * L_u * Deg^(-0.5);
% 
% K_u = inv((eye(N,N) + sigma * L_u_reg));
% K_u_inv = inv(K_u);

% save('u1_user_regLap_kernel.mat', 'K_u');
% save('u1_user_regLap_kernel_inv.mat', 'K_u_inv');



